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ABSTRACT 

We present a fully nonlinear hydro dynamic 'shallow water' model of the solar 
tachocline. The model consists of a global spherical shell of differentially rotat- 
ing fluid, which has a deformable top, thus allowing motions in radial directions 
along with latitude and longitude directions. When the system is perturbed, 
in the course of its nonlinear evolution it can generate unstable low-frequency 
shallow-water shear modes from the differential rotation, high-frequency gravity 
waves, and their interactions. Radiative and overshoot tachoclines are character- 
ized in this model by high and low effective gravity values respectively. Building 
a semi-implicit spectral scheme containing very low numerical diffusion, we per- 
form nonlinear evolution of shallow- water modes. Our first results show that, (i) 
high latitude jets or polar spin-up occurs due to nonlinear evolution of unstable 
hydrodynamic shallow-water disturbances and differential rotation, (ii) Reynolds 
stresses in the disturbances together with changing shell-thickness and meridional 
flow are responsible for the evolution of differential rotation, (iii) disturbance en- 
ergy primarily remains concentrated in the lowest longitudinal wavenumbers, (iv) 
an oscillation in energy between perturbed and unperturbed states occurs due to 
evolution of these modes in a nearly dissipation-free system, and (v) disturbances 
are geostrophic, but occasional nonadjustment in geostrophic balance can occur, 
particularly in the case of high effective gravity, leading to generation of gravity 
waves. We also find that a linearly stable differential rotation profile remains 
nonlinearly stable. 

1. Introduction 

The solar tachocline plays a fundamental role in governing the dynamics and evolution 
of the Sun's interior magnetic fields, manifestations of which are observed at the solar surface. 
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The tachocline is a very thin (~ 0.03-R) transition layer at the base of the convection zone 
through which the solar differential rotation changes from latitudinal in the convection zone 
to a uniform core rotation. This layer contains strong radial and weak latitudinal differential 
rotation. 

Due to existence of both radial and latitudinal shear in the tachocline, the study of its 
dynamics and physics is complex. Analysis of tachocline hydro- and magnetohydrodynamics 
have proceeded in two different ways: the first way treats the thin tac hocline explicitly as a 



full th ree dimensional multilayered spherical system (see, for example, IChan. Liao k Zhang 



( 120081 )). The other way is to treat the tachocline as a one-layer system. The multilayer 
system has been used to study the dynamics of both shears, radial and latitudinal, while 
the one-layer system has been applied to examine the dynamics and stability of only the 
latitudinal shear. Over the past decade, various models have been built for global HD/MHD 
of the tachocline latitudinal shear - 2D, quasi-3D shallow- water systems and 3D thin-shell 
primitive equation systems, under both the Boussinesq and hydrostatic approximations. 

The linear and nonlinear calculations of HD/MHD instability of tachocline latitudi- 
nal shear have been performed in great detail for a 2D (latitude-longitude) system. These 
studies led to the following major outco mes: (i) solar tachocline latitudinal shear is hydro- 



dynamically stable in pure 2D systems (IWatson I Il98ll ; iDziembowski k Kosovichev 1 11987 



Charbonneau. Dikpati k Gilman 111999c iGaraud 11200 ih : (ii) MHD shear instability with lon- 
gitudinal wavenumber m = 1 to 7 sets in and persists in the solar tachocline in the presence of 
a wide ra nge of toroidal field profiles (from broad to narrow bands), amplitudes and latitude- 
locations fiGilman k Fox1ll997l : Dikpati k Gilman Ill999[ iGilman k Dikpati 1120021 1: (iii) ax- 
isymmetric HD/MHD instability does not, or more precisely, cannot exist in 2D; (iv) as 
consequences of the no nlinear evolu tion of this instability, broad t oroidal fields open up into 



" clam -shell" patterns ( jCally 1 120011 ). whereas narrow bands tip ( jCally. Dikpati k Gilman 
20031 ). 



We now know that the 2D results are robust. All of them obtained fro m linear analysis 
have been found to appear in more complex quasi-3D 'shallow-water ' systems (IDikpati k Gilman 
200ll : IGilman k Dikpati I I2OO2I bikpati. Gilman k Rempel I l2003h . in which the third di- 
mension is generally included in a restricted way, namely the upper boundary of the shallow 
fluid-shell is allowed to vary with latitude, longitude and time so that the radial motions can 
occur. Intuitively it is expected that the inclusion of the third dimension would produce some 
other unstable modes that are purely characteristic of 3D because, in addition to the kinetic 
and magnetic energy reservoirs due to the presence of differential rotation and magnetic 
field, 3D systems allow energy extraction from the potential energy reservoir, particularly 
when the effective gravity (G) of the system is low. G is a measure of the subadiabatic 
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stratification of the fluid layer; low G can allow more deformation of the fluid-layer's top 
surface. In linear studies, it was indeed found that a hydrodynamically stable 2D tachocline 
can be unstable to quasi-3D per turbations, due to the p resence of this additional energy 
source in a shallow- water model (IDikpati & Gilman Il200lf ). 



Many idealized shallow-water models have been built over the past century, since the 
first, single-layer model by Hough (1898). For studying the mid-latitude oceanic tides and 
atmos pheric zonal flows and meridional circulation , mult i-layer shallow-water models are 
used ( iHuang fc Pedlosky Ill998l ; iHurlburt fc Hogan 1 120081 ) . Shallow water models are a lso 
being used for operational forecasts of El Nino and La Nina (IPhilander &: Fedorov 1120031 ) . 



Shallow water theory was applied to the global solar tachocline by Gilman (2000), who 
first formulated an MHD analog of a shallow-water system. Since then, all global studies 
of HD/MHD shallow- water instabilities in the solar tachocline have been performed in the 
linear regime. Nonetheless, even within the scope of a linear model, the upward bulging 
of a magnetized shallow fluid layer created at certain longitudes, by growing modes with 
low-longitudi nal wavenumbers, can prov ide a plausible explanation of the Sun's "active lon- 



gitudes" (see IDikpati fc Gilman I ( 120051 ) for details). This linear model was able to simulate 



active longitudes up to 15 Carrington Rotations, after which nonlinear interactions among 
the unstable modes start playing important roles. In order to be able to simulate active 
longitudes for one full solar cycle and several solar cycles, it is necessary to study the the 
nonlinear evolution of the aforementioned instability in the tachocline fluid layer. 

Over the past decade, several 3D thin-shell and fu ll 3D HP and MHD models for linear 



analysis of tachocline shear instab i lities have been built (ICa lly 



20031 : lArlt. Sule & Riidiger 1 120051 : iGilman. Dikpati & Miesch 



2001 



Zhang. Liao fc Schubert 



20071 ). However, unlike many 



studies of nonlinea r evolution of sh a llow-w ater instabilities in the oceanic and atmospheric 



context (see, e.g. iPoulin fc Flierl I (120031 )). there exist only a few studies of the nonlin- 



ear evolution of 3D MHD instabilities in the solar contex t, using 3D thin-shell models 
( IMiesch. Gilman fc Dikpati 1 120071 : iHollerbach fc Cally 1120091 ). Although in a shallow- water 
model the radial dimension of the fluid layer is included in a simplified way by allowing the 
layer's thickness to vary with latitude, longitu de and time, it is no t as re stricted compared 
to a full 3D MHD system as it might appear. iRempel &: Dikpati I ( 120031 ) have shown that 
an MHD shallow-water approach is a first-order approximation of a full MHD equilibrium 
in which the magnetic curvature stress is balanced by a latitudinal pressure gradient for a 
strongly subadiabatic stratification and by modification in the shear flow for a weakly suba- 
diabatic stratification. The height deformation of the top surface layer in the shallow-water 
approach corresponds to the deformation of the surfaces of constant entropy required for 
latitudinal force balance in the full MHD approach. Thus using a shallow water system for 
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the Sun, specifically for the solar tachocline, can help reduce the complexity of a full 3D 
MHD system without much affecting the basic physics behind the global dynamics of the 
fluid layer. 

Putting aside the added complexity due to the presence of magnetic fields, we first focus 
on the nonlinear evolution of shear instability in an unmagnetized solar tachocline and find 
out how the tachocline latitudinal differential rotation profile can evolve with time and get 
modified due to the influence of this instability. In pure 2D calculations, the jets were seen 
to form around the locations of the singular points arising from the coincidence between the 
Doppler-shifted phase-speed of the unstable modes and the she ar-speed of the system, or 



the A lfvenic singular points in the case of magnetized tachocline ( ICally. Dikpati fc Gilman 



20041 ) . Along with the exploration of jet-formation on unperturbed differential rotation 
profiles due to the nonlinear evolution of hydrodynamic shear instability, another major aim 
of this paper is to study the nonlinear evolution of the bulging and depression of the fluid 
layer in conjunction with the shear flow pattern, or in other words, the exchanges between 
the kinetic and potential energies of the system. 

We will present the mathematical formulation including equations and solution tech- 
niques in §2 and in the Appendix, and the implementation of spectral viscosity in §3. In the 
description of results in §4 we will present in the first subsection how the energy spectrum 
develops, i.e. how the total energy gets distributed in different longitudinal wavenumbers 
and how they evolve with time. The subsequent subsections of §4 will focus on the evolu- 
tion of the Sun's longitude-averaged differential rotation and the height profile, the energy 
exchanges between kinetic and potential energy reservoirs, and the evolution of flow and 
deformation of the fluid layer on the latitude-longitude surface. We will summarize our 
findings in §5. 



2. Mathematical Formulation 



2.1. Governing equations 



Hydrodynamic shallow-water equations were built for studying global atmos pheric and 



ocean ic dynamics; they exist in the literatur e in va rious forms, for example see iPedlosky 



(119871 ) . Following that, iDikpati fc Gilman I (120011 ) described the hydrodynamic shallow- 
water equations, along with the usual assumptions and approximations in their §2.1 and 
§2.2 and solved them in the linear regime. Briefly, in a shallow water model the velocity (V) 
of the fluid can be defined as V = uX + v<f> + wr, where u and v are the horizontal velocity 
components in longitude (A) and latitude ((f)) respectively, and w is the radial velocity. In the 



shallow water approximation, u and v are independent of depth, and w is a linear function 
of depth. For a shallow fluid layer with a rigid bottom and a deformable top, the maximum 
vertical velocity will be at the top. 



Due to the thinness of the tachocline compared to the solar radius, the divergence of 
radius and the density variation within the shallow fluid layer are ignored in the momentum 
and mass-continuity equations. All time scales of the system are considered much longer 
than the acoustic time scale. With these assumptions, the hydrodynamic shallow water 
equations for a sin gle-layer tachocline can be derived from the mass continuity and momen- 
tum equations (see iDikpati fc Gilman I ( 1200 ll ) for detailed derivation). In order to study the 
nonlinear evolution of the hydrodynamic shallo w- water instability i n the solar tachocline, 
we start from the equations (4), (5) and (8) of IDikpati &: Gilman I ( 120011 ) and rewrite the 
nonlinear dimensional hydrodynamic shallow-water equations in the rotating frame as, 
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in which u and v are the horizontal velocity components in longitude (A) and latitude (0), 
and h is the deformable top surface. 

We make equations (1-3) dimensionless by choosing the radius of the shell (r ) as the 
length-scale and the inverse of the core rotation rate (w c _1 ) as the time-scale of the system; 
thus the thickness of the fluid layer representing either the radiative or the overshoot part of 
the solar tachocline should be between 0.01 and 0.05 dimensionless units, and one year time 
will correspond to about 100 dimensionless units. The nondimensional equations are 
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in which G = gHjr\uj\ is the effective gravity. With typical values H, the pressure scale 
height, the value of G is high (> 10) in the radiative tachocline and low (< 1) in the overshoot 
part of the tachocline. 

Due to the pole problem, it is difficult to solve the advective form of shallow-water 
equations in spherical polar coordinates. It is a fundamental problem for solving numerically 
the partial differential equations on the sphere, particularly using finite difference methods - 
certain terms are likely to be unbounded in the neighborhood of the poles. Spectral methods 
are popular for their high accuracy and their ability to avoid the pole problem. There exists 
a vast literature on h ow to avoid the pole problem. We follow here a technique proposed by 
Swarztrauber I (119961 ) and so we define the vorticity (() and divergence (5) as 
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and rewrite the equations (4-6) as follows: 
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where,/ = 2u c sin is the Coriolis parameter. 



We apply the spectral transform method ( ISwarztrauber Ill996l ) to solve the Equations 
(8-10), by expanding the vector field (i.e. the flow) and the scalar field (height) respectively 
in terms of vector and scalar harmonics. Thus the pole problem can be handled efficiently 
without raising the order of the equations. The appendix describes in detail the technique 
used to solve the Equations (8-10). 



2.2. Description of Initial State, Numerical Algorithm and Spectral Viscosity 



We express the solar tachocline latitudinal differential rotation in angular measure, oj, 
so that u — oj cos (p. From the knowledge of helioseismic observations the mathematical form 



-7- 



of lo has been adopted as 



u = s — s 2 /i 2 — s 4 /i 4 — ui c , (11) 



in which w c represents the solid rotation of the core, and Sq, s 2 and S4 are numerical coeffi- 
cients. (S2 + S4)/so ~ 0.3 represents solar latitudinal differentia l rotation near the surface, 

wher eas in the radiative part of the tachocline, (S2+S4) / so < 0.17 (jcharbonneau. Dikpati fc Gilman" 



19991 ). Detailed parameter surveys in the linear regime (IDikpati &: Gilman 1 120011 ) have re- 
vealed that oj becomes unstable to perturbations with longitudinal wavenumbers m = 1 and 
m = 2 for a wide range of s 2 and s 4 values. In order to focus our study on how this insta- 
bility evolves in the cases of high and low effective gravity (G), we pick an u with s 2 = 0.15, 
54=0.09 and So = 1.051 so that the pole-to-equator differential rotation is 24%. Linear 
studies indicate that this differential rotation is hydrodynamically unstable for a wide range 
of G. Therefore we can demonstrate the features of the time-evolution of this differential 



rotat i on for both high G and low G case s, in particular G = 10 and 0.5 (IRempel fc Dikpati 



2003c IDikpati. Gilman &: Rempel 1 120031 ) which are, respectively, characteristic of radiative 



and overshoot tachoclines. 

We also study a few other cases with 30% and 15% pole-to-equator differential rotation. 
In all cases we start with an initial reference state profile that is a combination of unperturbed 
profile and a perturbation with low longitudinal wavenumbers (mostly m = 1 mode), with 
a specified amplitude with respect to the unperturbed differential rotation. 

In all simulations we start from known solutions to the linearized instability equations, 
rather than from perturbations of arbitrary form, which has the effect of initially filtering 
gravity waves out of the calculations. But once nonlinearities become important, gravity 
wave modes can and frequently are excited and subsequently participate in the overall dy- 
namics. 



2.3. Spectral Viscosity 

Spectral methods, despite their accuracy and absence of pole problems, are prone to 
generating spurious oscillations in physical space and hence a reduction in global accuracy. 
In the case of nonlinear evolution equations in an idealized, inviscid model like the one we are 
presently dealing with, spectral methods generate energy cascade from lower to higher wave 
number modes, which often creates a spectrum with an erroneously high amplitude tail at 
high longitudinal wavenumbers - the so-called spectral blocking problem, which could even- 
tually lead to a poor energy-conservation during the nonlinear evolution of the system. The 
traditional approach is to use an artificial horizontal diffusion to deal with these difficulties 
of spectral methods. 
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Fig. 1. — Spherical spectral viscosity operator as function of wavenumber I in Legendre 
polynomial P™. 



We apply the Spectral Viscosity (SV), as prescribed by iGelb &: Gleeson I (j200l[ ). which 
is a spherical spectral viscosity operator, able to conserve invariants of the system at no 
additional computing cost and to retain the high wave number information in the system. 



The dynamical variables a ™, a™, b™, b™, c ™ and q m are filtered at each time step using 
the formula 
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in which M is the truncation limit of n. The e and l c have been taken to be e 
l c = 2M 3//4 . The pattern of af v as function of I is shown in Figure 1. 



m and 



-9 - 



In a viscous model for the astrophysical fluid simulations, the use of artificial horizontal 
diff usion will not be needed, as it would always include the physical diffusion (for example, 



sec 



Fan et al. I (119991 )). The comparison between an inviscid model with artificial diffusion 
and a viscous model without artificial diffusion can be done by using the similar amount of 
artificial and physical diffusion in respective models. In order to derive an equivalent physical 
diffusion from the artificial diffusion used in an inviscid model, a plausible experiment can 
be set up as follows. By switching off all energy sources, such as kinetic energy from shear 
flow and the potential energy from effective gravity, just the time-rate of the spread of a 
2D Gaussian function in latitude-longitude direction can be studied in the presence of an 
artificial diffusion. The horizontal diffusion, estimated from the rate of spread of 2D Gaussian 
function, can be used as the amount of physical diffusion in a viscous model without any 
artificial diffusion. The nonlinear evolution of a shallow-water instability, specifically the 
spectral blocking issue and the quality of energy conservation from such a viscous model can 
be compared with that from an inviscid model with artificial diffusion. A forthcoming paper 
will address this comparison quantitatively. 

However, one of the limitations of a single-layer shallow-water model, in both viscous 
and inviscid cases, is that the horizontal velocities are independent of depth - only the 
vertical velocity varies with depth. Thus there are no viscous terms that involve radial 
second derivatives of horizontal velocities in the horizontal components of the momentum 
equations. In reality, a full 3D tachocline model will be a viscous model. So a multilayer 
shallow-water model will be needed to connect the layers through some diffusion in the 
radial direction, as has been done in the case of a full 3D, multilayer tachocline model by 



Chan. Liao fc Zhang I (120081 ). Alternatively, the effect of radial diffusion can be explored by 



using a Newton's cooling type for malism in which the radial diffusi on terms can be converted 



into drag terms (for example, see iDikpati. Cally fc Gilman I ( 120041 ) ). In future, development 



of such a model will also help investigate the effect of forcing the differential rotation at the 
top of the layer. 



2.4. Adjusting the time step 

In order to increase the efficiency of the numerical code, we use an adjustable time 
step which can be changed (up or down) during the the simulation. If the time step, At, is 
taken large, it may result in poor energy conservation, whereas if it is taken small, it will 
lead to longer runs that will take many steps and more time to compute. In the course of 
a computation, if the system is going through an interval near the peak of its potential or 
kinetic energy, then a somewhat lower time step will ensure better accuracy. At other times 
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of the evolution, a larger time step can be taken, because that will speed up the computation. 
Therefore we have implemented a variable time step that is adjustable at each iteration. The 
scheme for adjusting time step is described below. 

We define e as the fractional update of a variable- the ratio of the magnitude of the 
update of that variable and the magnitude of the variable at that time step. When evaluating 
e we exclude those variables whose absolute magnitude is less than 10~ 4 , because the effect 
of those variables on the physical dynamics at that point in time is very small compared to 
the normalized variables. Then we evaluate the maximum value of e over all the updates at 
that time step. If the maximum in e turns out to be less than 0.005 then we scale up the next 
time step by a factor of 1.02 of current time step. On the other hand, if the maximum in e 
is greater than 0.005 then we scale down the time step by 0.8 and reevaluate the evolution 
of the variables at that time step. In this way we limit the amount of maximum update of 
a variable, and hence the drastic change in the system at a given time step, to achieve good 
energy conservation over long periods of time. 

Values of scale-down and scale-up in the time-step are chosen to be 20% and 2% re- 
spectively. These are somwhat arbitrary, but we have found for our present calculations 
the amount of scale-down in the time-step needs to be larger than that of scale-up. This is 
because accuracy in the conservation of the energy deteriorates when the system is evolving 
rapidly. 

3. Results 

In order to present results and physical interpretations in a well-organized fashion, we 
divide this section into four subsections, which will successively present the convergence test 
related to the selection of longitudinal wavenumbers, a test for total energy conservation and 
exchanges among different energy reservoirs, nonlinear evolution of longitude-averaged flow 
and height profiles, and nonlinear evolution of flow accompanied by deformations of the fluid 
shell in latitude and longitude. 

3.1. Convergence test for truncation in m 

For the technique described in detail in the appendix, the computation speed largely 
depends on the number of Z's and m's we select in the basis functions. The total energy is 



defined in this system as 
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in which u, v and h represent longitudinal velocity, latitudinal velocity and height deforma- 
tion. Since the modes with longitudinal wavenumbers m — 1 and 2 are the only modes that 
are linearly unstable in the solar tachocline, we construct our initial configuration of the flow 
and height-deformation by using the modes with m equal to 2 at the most. 

We calculate the energy distribution at each longitudinal wavenumber m, using the 
expression (12), for two G values, namely G = 10 and G = 0.5, and plot them in Figures 
2(a) and (b) as functions of m. In both cases we have taken a pole-to-equator differential 
rotation amplitude of 24% with s 2 = 0.15, s 4 = 0.09 in the expression (13). The total energy 
of the system for different wavenumbers m = 1 through 10, plotted at the time t — 0, t — 50 
and t — 90 (respectively by black, red and blue diamonds) show that the energy contribution 
falls off rapidly with m. At t — the total energy of the system is contained in m = and 
m — 1, because that is how we constructed the reference state, which suggests that m may 
be truncated at 6 or 7 with very little loss of accuracy; this speeds up our computation 
enormously. 



-12- 




m 



Fig. 2. — Spectrum of total energy as function of longitudinal wavenumber m at t = (black 
diamonds), 50 (red diamonds) and 90 (blue diamonds), for (a) G = 10 and (b) G = 0.5. 
Note that dimensional time for t = 40 and 90 are respectively 4.8 and 10.8 months. 



3.2. Evolution of energetics 



Because the fluid-shell system with latitudinal differential rotation has its own inherent 
instability for low-order longitudinal wavenumber modes, in the course of the nonlinear 
evolution of such a system its unstable modes grow by taking energy from the unperturbed 
system and thus modifying the system. At a certain point the reference-state system gets 
modified to the extent that it takes back energy from the perturbations. The cascading 
of energy to higher longitudinal modes is not significant in this system for m below 6 or 7; 
instead a stable oscillation of energy between the reference-state and perturbations develops. 
Figures 3(a) shows such an oscillation between the reference-state (m = 0) kinetic energy 
(K) and the perturbation (m > 0) kinetic energy (K'), for a perturbation of about 40% in 
the differential rotation (angular measure, i.e. u = ucos8), for the case of G = 0.5. Similar 
oscillations occur in G = 10 case. 

The same period and phase of oscillation occurs between the reference state (P) and 
perturbation (P') potential energies in the system, but with substantially smaller amplitude 
than for the kinetic energy. We see that during this oscillation the total energy of the 
system remains virtually constant, which shows that the numerical viscosity in the system 
is extremely small. In fact, the total kinetic and total potential energies are separately 
constant, so that the energy exchanges are strictly between reference state and perturbation 
energies of the same type. 
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Fig. 3. — Frame (a) shows the evolution of reference state kinetic energy, K (i.e., KE m= o, 
plotted in red solid curve), potential energy, P (i.e., PE m=0 , blue solid curve), perturbation 
kinetic energy, K' (i.e., KE m>0 , red dashed curve), and potential energy, P' (i.e., PE m>0 , 
blue dashed curve). The energy exchanges between K and K', and between P and P' shows 
an oscillatory pattern, the frequency of which depends on the amplitude of the perturbation, 
while total energy, plotted in solid black curve, is conserved. Frame (b) shows the oscilla- 
tion period of energy exchange between the reference state and perturbation as function of 
perturbation energy measured as a fraction of the reference state energy of the system. 
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From equation (56) of Dikpati & Gilman (2001) we can infer that this oscillation is 
driven by the time dependent Reynolds stress of the perturbations. Initially this stress 
extracts kinetic energy from the differential rotation, which is modified in such a way as 
to reduce the energy available for continued growth in the disturbance. But the system 
then 'overshoots' an equilibrium state, resulting at later phases in the oscillation in the 
perturbations giving back kinetic energy to the differential rotation. Then the perturbations 
are able to grow again, and the oscillation repeats. In §3.4 we will show how the disturbance 
structure changes during the oscillation. 

Figure 3(b) shows how the oscillation period depends on the initial perturbation am- 
plitude. The relationship is obviously not linear. The process that causes this behavior is 
also nonlinear. Physically what is happening is that as the initial perturbation amplitude is 
increased, energy is extracted from the differential rotation at a faster rate. This modifies 
the differential rotation profile faster to a form that is less unstable to the perturbation and 
subsequently to a form that extracts energy from the perturbation. This makes the period 
shorter. There is probably a limit to how short the period can be, as evidenced by the curve 
in Figure 3b approaching an asymptote in period that is greater than zero. 

Conversely for very low initial amplitude the period gets longer, but probably not to 
infinity because the initial growth is exponential (and the rate is independent of the small 
amplitude). It simply takes longer for the growing disturbance to reach a point in time at 
which the feedback from the modified differential rotation is strong enough to limit further 
growth. 

3.3. Nonlinear evolution of differential rotation formation of high-latitude 

jets 

We have seen above that the kinetic and potential energies of the nonlinear shallow 
water system go through a well defined oscillation. How is this manifested in the changes in 
the differential rotation profile? Figures 4a,b show its behavior for the same time period, for 
both high G (frame a) and low G (frame b), for a 40% initial disturbance amplitude. Here is 
plotted the angular velocity in the shell for five later times, compared with the initial profile 
(black curve). All angular velocities are taken relative to the core rotation rate. 
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Fig. 4. — Frame (a) shows, for a high effective gravity (G = 10) snapshots of longitude- 
averaged differential rotation at several selected time during the course of its evolution. 
Frame (b) shows the same, but for a low effective gravity (G = 0.5). 
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We see in Figure 4 that the changes in angular velocity with time are much larger in high 
than in low latitudes. This is because it is angular momentum that is being transported in 
latitude by the Reynolds stresses of the perturbations, and the moment of inertia and volume 
of high latitudes is much smaller than that of low latitudes. So depositing a given amount 
of angular momentum in high latitudes causes the angular velocity to rise much more there 
than it falls in low altitudes. It follows that if this instability is active in the solar tachocline, 
the observed high latitude angular velocity changes would be much larger than those in low 
latitudes. To date, helioseismic measurements have not been able to tell us whether this is 
so, because high latitudes are very hard to observe. 

For both high and low G, the peak high latitude angular velocity is reached within a 
little over one month after the initial disturbance is introduced. Thereafter it fluctuates, but 
about a distinctly higher polar angular velocity than in the initial state. The polar angular 
velocity is significantly higher for low G than high. We reason below that this is because with 
low G it is easier for the mass distribution in the shell to be rearranged by meridional flow. 
The detailed time history of the amplitude of the polar angular velocity does not precisely 
match the rather smooth oscillation seen in Figure 3. This is not surprising, because the 
time variation seen in Figure 3 is of global integrals of the energy, while the polar rotation 
changes are for a very small part of the volume of the global shell. 

We can illustrate and explain the physical processes involved in the evolution of the 
differential rotation by using equations (60)-(62) of Dikpati & Gilman (2001), which we 
repeat here, rearranged so that only the time derivatives are on the left hand sides. 

^ = ^-^(co(l-/, 2 )), (13) 
^ = F v - 2 u W 2f x + Gy/T^tf^ (14) 

^ = f h - A (a + m wr 3 ^)) , (15) 

The forcing functions F u , F v and F h are defined as, 



F u = -V'^u'y/l-fj?), 



Fk = — (wyry). 
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Here the variables are H2, V2 and h 2 , which represent respectively the departures of linear 
rotational velocity, meridional flow, and layer thickness from their initial values, uo is the 
angular velocity of the core, and /1 is cosine of the latitude. From above definitions, one can 
see that F u , F v and Fh are quadratic functions that contain the covariances or stresses from 
the perturbation variables u, v, h. Equations (13) and (14) are, respectively, the equations of 
motion for differential rotation and meridional flow, and Equation (15) is the mass continuity 
equation, all for the longitudinal wavenumber m = variables. 

Figure 5 displays profiles of layer thickness, meridional flow and angular velocity at 
the same times as shown in Figure 4b. Thickness (solid purple curve) and latitudinal flow 
(solid green curve) are shown in frames (a) through (e), along with the initial thickness 
(dashed purple curve) for reference. Note that initial meridional flow (dashed green curve) is 
zero. The differential rotation angular velocity (thick black curve) is shown in frames (f)-(j), 
together with the initial differential rotation (thin curve). 

If we compare the thickness and latitudinal flow profiles as time advances, we can see 
clearly that initially the thickness shrinks near the poles but later rises again as the oscillation 
progresses. Consistently, the meridional flow is away from the poles while the thickness is 
declining there, and toward the poles while the polar thickness is rebounding. The same 
relationship between meridional flow and thickness, dictated by mass conservation, can be 
seen at all other latitudes. When the polar thickness is declining, the polar angular velocity 
is increasing. 

The physical processes that govern these correlated changes among thic kness, meridional 



flow a nd differential rotation are contained in Equations (13)-(15) above (see lDikpati fc Gilman 



( 1200 11 )). In brief, what happens during the nonlinear oscillation is that, through F u (Equa- 
tion 13) the perturbation Reynolds stress transports angular momentum toward the poles to 
cause a 'spin-up', or higher linear and angular velocity there. In the meridional flow Equa- 
tion (14), this 'spin up' leads to a negative (equatorward) Coriolis force that pushes fluid 
away from the pole toward the equator. This mass flow then produces a positive latitudinal 
pressure gradient force, which is expressed as a thickness gradient in a shallow-water model. 
This force starts to counterbalance the Coriolis force. 
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Fig. 5. — Frames (a-e) show, for G = 0.5, five snapshots of longitude-averaged height defor- 
mation (solid purple curves) for the same selected times as in Figure 4(b), and corresponding 
latitudinal flow (solid green curves) have been superimposed on that. Dashed purple and 
green curves in each frame respectively present height deformation and latitudinal flow at 
t = 0. Solid black curves in frames (f-j) show longitude-averaged differential rotation, for 
the same selected times, with a superimposed dashed curve for t = 0. The right column is 
created by separating five curves of Figure 4b. 
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This rebalancing of forces will reduce the equatorward meridional flow, but in our nearly 
dissipation free system, it overshoots, causing the meridional flow to reverse and rebuild the 
polar thickness. At the same time, the changed differential rotation profile is less unstable to 
the perturbations, or even stable to them, so the Reynolds stresses get much weaker, and the 
polar latitudes spin down again. Throughout this sequence of changes, the mass continuity 
Equation (15) ensures that the changes in thickness are consistent with the changes in 
meridional flow. 



3.4. Evolution of disturbance planforms 

We can better understand the nonlinear dynamics of the shallow water system by ex- 
amining the planforms (longitude-latitude) of the flow and thickness as they evolve. Figure 
6 displays these quantities for the same times as in Figures 4b and 5 (low G, 0.5). The left 
hand column displays the total flow and total thickness fields (blue is thinest, red thickest); 
the right hand column shows all of the m > fields, after subtracting out the axisymmetric 
parts of the differential rotation, meridional flow and axisymmetric thickness. 

From both columns we can see that the m — 1 mode dominates in the flow and thickness 
patterns, leading to a 'meandering' of the E-W flow. As time progresses, more structure 
and somewhat finer spatial scales develop, as we should expect from the behavior of the 
energy spectrum shown in Figure 2. We can also see from both columns that the flow is 
predominantly geostrophic, that is in the North there is counterclockwise circulation around 
the blue (thinest) areas, clockwise around the red (thickest) areas; the reverse is true in the 
South. Since the pressure in the shallow water system is hydrostatic, it is proportional to the 
thickness, so the counterclockwise flow in the North is around low pressure, the clockwise 
flow around high pressure. 

From the right hand column we can see that in the first 1-2 months the disturbances 
are transporting angular momentum toward the poles to spin them up, as seen in Figure 
4b. The velocity and thickness patterns are tilted with respect to the E-W direction in such 
a way that the fluid elements moving toward the east in the disturbance are also moving 
toward the poles, while those moving toward the west (relative to the rotating reference 
frame) are moving toward the equator. This means there is a correlation between the two 
velocity components that is the Reynolds stress, which is carrying out the poleward angular 
momentum transport. 
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Fig. 6. — Frames (a-f) show, for G = 0.5, planforms of flow (arrow vector) and height (color 
map) for five selected times as in Figure 5. Frames (g-1) show the perturbation part of flow 
and height. Red represents swelling and blue depression. 
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Fig. 7. — Same as in Figure 6, but in Mollweide projection. 
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By 2.1 months, these tilts are virtually gone, indicating angular momentum is no longer 
being transported toward the poles. This is consistent with the profiles of differential rotation 
shown in Figure 4, in which the maximum amount of polar spin- up is reached between 1.3 and 
2.1 months in these selected cases. By six months, the perturbation patterns have become 
more complex, and it is difficult to discern which way the angular momentum transport 
is going. But from Figure 3 there is still a pronounced amplitude oscillation between the 
perturbations and the differential rotation and meridional flow. In this nearly dissipation- 
free system there does not appear to be an absolutely steady nonlinear state. Instead there 
is a dynamical equilibrium with continuous evolution of all flow patterns within a nonlinear 
regime in which the high latitude rotation rate, while fluctuating, is almost always faster 
than in the initial state. 

In Figure 7, we repeat the twelve frames (Figure 6a-l) in Mollweide projection. This 
Figure captures a more realistic look of the global disturbance patterns, particularly at and 
near the poles. In latitude-longitude plots, the single point at each pole is spread out over 
360° longitude. Figure 7 reveals all the features of the disturbances seen in Figure 6, but 
also reveals more clearly that the flow spreads over both hemispheres crossing the equator. 

3.5. Is a linearly stable profile nonlinearly unstable? 

To illustrate further the range of nonlinear behavior of the shallow water system, we 
have performed an additional numerical experiment in which we have chosen a combination 
of differential rotation (s = 0.15) and effective gravity (G = 10) that we know from linear 
theory is stable to perturbations, and initially perturbed the system with a disturbance 
calculated for an unstable case. The amplitude of perturbation used is 20% of the reference 
state differential rotation. Shown in Figure 8 is the resulting differential rotation evolution 
sampled at the same times as in previous cases. Here we see that in the first few days of 
the integration, the disturbance starts to spin up the poles (red curve) because the Reynolds 
stress in it transports angular momentum toward the poles. But this trend is very short-lived, 
as the stable differential rotation modifies the disturbance, changing the sign of its Reynolds 
stress by changing the sense of tilt of the perturbation velocities with respect to the E-W 
direction. Consequently for the next two months of the integration, the poles actually spin 
down to lower rotation than they had initially, (light green curve) while the kinetic energy 
of the differential rotation is increased at the expense of the energy in the disturbance. But 
this, too, is a transient, and by seven months the rotation at high latitudes has returned to 
very close to its initial profile (dark blue curve). 

The details of how the disturbance evolves with time is shown in Figure 9, which gives 
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Fig. 8. — Response of a linearly stable solar differential rotation (with s = 0.15 and for G = 
10) to perturbation. Five snapshots present the evolution of longitude-averaged differential 
rotation for five selected times as in Figure 5. 



longitude-latitude planforms of velocity and shell thickness. The three frames are for, re- 
spectively, t = 0, 10.6 and 49.8 (0, 1.3 months, and 6 months). We see that in this succession 
of frames the tilt of the velocity vectors with respect to a latitude circle changes orientation 
twice, as the associated Reynolds stress changes sign twice, signifying first poleward, then 
equatorward, then poleward momentum transport again. This is what causes the poles to 
spin up, then spin down, then up again. 
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Fig. 9. — Frames (a-c) show, for s = 0.15 and G = 10, planforms of flow (arrow vector) and 
height (color map) in Mollweide projection, for three selected times, t = 0, 1.3 months and 
6 months. 
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Since linear studies ( IDikpati &: Gilman 1120011 ) indicate that the low differential ampli- 
tudes (s < 0.17) are stable in the radiative tachocline with high G, we considered a case 
with s — 0.15 and G = 10 in this subsection. In linear studies shallow-water gravity modes 
were extracted out from the calculation and only shear modes were explored. However, we 
have no such freedom in the case of nonlinear evolution of a shear profile in the radiative 
tachocline fluid shell with high G. Consequently we see the appearance of gravity modes in 
the evolution discussed in Figures 8 and 9. The generation and evolution of gravity modes is 
a vast subject of research itself. A forthcoming paper will address the generation of shallow- 
water gravity modes due to nonadjustment of geostrophic balance when the effective gravity 
of the system is high. 



4. Concluding Comments 

We have presented above the first fully nonlinear 'shallow water' numerical model to be 
applied to the Sun. It is intended for use in understanding the global dynamics of solar and 
stellar tachoclines. It provides a beginning and template for later models that will include 
global MHD effects; solar tachoclines most likely contain strong global scale magnetic fields, 
particularly toroidal fields. 

This spectral model needs only extremely low numerical diffusion to keep the solutions 
well behaved, so it is particularly well suited to simulation of strongly nonlinear time depen- 
dent hydrodynamics and MHD of the Sun. The equations conserve energy with extremely 
high accuracy (~ 0.001%), so that very long simulations are possible without significant 
spurious build-up of kinetic or potential energy at small scales. Substantial increases in both 
computational efficiency and accuracy are obtained by use of a mechanism for adjusting the 
time step during a simulation. 

Since the model does include the radial dimension in a simple way, by means of a 
deforming top surface, it contains gravity wave modes as well as shearing instabilities of the 
differential rotation. It also allows for a simplified form of meridional flow that is coupled to 
and causes the changes in the thickness of the spherical shell. Thus it is possible with this 
model to study rather complex fluid dynamics, including interplay between the relatively 
low frequency global shearing flows and instabilities, and the higher frequency gravity waves 
made possible by the deforming top boundary. 

The low dissipation in the model is of particular value for simulating the global dynamics 
of the 'radiative' tachocline that is found below the overshooting layer of the convection zone 
above. But it is equally valuable for simulations applied to the 'overshoot' tachocline, whose 
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differential rotation is maintained by downward diffusion from the convection zone. 

Our first results from this nonlinear model focus primarily on the nonlinear interac- 
tions between unstable global perturbations previously found from linear analyses and the 
differential rotation that gives rise to them. We show that the nonlinear dynamics result 
in a 'spin-up' of polar regions due to the Reynolds stress associated with the growing dis- 
turbances, and also substantial nonlinear oscillations between disturbance and differential 
rotation kinetic energy. The oscillation occurs because as the disturbances extract energy 
and transport momentum to high latitudes, the differential rotation profile changes in such 
a way as to reduce and even shut off the instability. At least temporarily there can occur a 
reversal in direction of angular momentum transport, which rebuilds energy in the differen- 
tial rotation. The simulations go through a sequence of these oscillations without settling to 
a steady state. This is because the system is virtually dissipation-free due to the implemen- 
tation of numerical scheme with very little spectral viscosity. Throughout this evolution, the 
longitude-latitude planform of the disturbances on the differential rotation, as well as the 
total flow, undergoes considerable evolution as well as propagation in longitude. 

The changes in angular velocity with time in polar latitudes are much larger than 
the compensating changes in low latitudes, because the moment of inertia and volume of 
polar latitudes is much smaller than in low and mid-latitudes. The period of the nonlinear 
oscillation is inversely proportional to the amplitude of the initial perturbations, apparently 
with finite asymptotic periods for both very large and very small initial disturbance energies. 
Throughout these oscillations the disturbance kinetic energy remains confined primarily to 
the lowest longitudinal wave numbers, m — 1,2, or the largest global scales. 

The simulations contain gravity waves of significant amplitude particularly when the 
effective gravity (G) of the system is high; their role will be addressed in a later paper. 
The solutions also contain kinetic helicity due to the correlation between the local radial 
component of vorticity and the radial motion associated with changes in the thickness of 
the shell. Kinetic helicity is of particular importance for the solar dynamo, and will also be 
studied in more detail in a later paper. In the Sun the differential rotation of the tachocline 
is imposed from the convection zone above. We will generalize our model to include this 
forcing in later studies. This will add still another time scale to the problem, namely the 
time it takes for the top forcing to restore the differential rotation. 

All of the additional studies just described are for the hydrodynamic model we have 
already built. The major step beyond that system is to include global MHD processes, 
which will without doubt substantially modify the hydrodynamic processes we have studied 
or will study soon. A still later step will be to couple this model to a global 3D flux 
transport dynamo model for the Sun, that will be capable of addressing the question of 
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the origin of active longitudes (IDikpati fc Gilman 1120051 ) as well as how such longitudinally 
dependent features and processes contribute the origin and evolution of solar cycles. The 
coupling would involve flow generated in the tachocline reaching into the convection zone 
above and generating 3D MHD induction effects. Thus there is a need to develop a viscous 
shallow-water model in order to be able to couple the tachocline instabilities with a global, 
flux-transport dynamo model in the solar convective envelope. In addition, it is necessary 
for the magnetic fields in the tachocline to connect with that in the convection zone above, 
continuously through the interface between the tachocline and the convection zone. To 
implement such a connection among magnetic fields at different layers, through the interfaces 
among the layers, we would follow the way as described in the figure 1 and equation (22) of 
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5. Appendix: Solution Technique 



In the case of primitive equations such as equations (8)-(10), the vector field (flow) 
and the scalar field (height) need to be expanded respectively in terms of vector and scalar 
spherical harmonics, in order to handle the pole problem. Pure scalar harmonics can also be 
used by converting all the vector fields into scalar fields, but at the expense of raising the 
order of the equations. 



Following the solution method of lSwarztrauber I (119961 ) (see method 2 there), the spectral 
decomposition of scalar h can be expressed as 



h(<f>,\) 
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and P\ m is the Associated Legendre polynomial of order (l,m). The * in the m summation 
implies a factor of | is multiplied when m equals 0. The coefficients aj" 1 and a' m can be 
obtained from the inverse transform as 

a 1 ™ = / d(f)P lm ((f)) cos0 / dX cosmXh ((f), X) , (A3) 
J-l Jo 

a[ m = [ 2 d(f)P lm ((f)) cos (f) [ dX sin mXh ((f), X) . (AA) 
J -I Jo 

In order to express the vector fields (u and v), we define the vector spherical harmonic 
components as 

d P 1 

y/l (l + Wm (0) = ^ = 2 K / + - m + l ) p Km-i) - Pl(m+1)} , (A5) 

y/l (l + l)W lm ((f)) = ^ = ±—L {P(,-i) (m+ i) + (/ + m)(Z + m - l)P(n) (m -i)} . (A6) 

cos<p 2 

The unit vectors can be constructed as 

iWi m \ - x _ ( —Vi m \ irn \ 



B lm =( "' m e*™* , C lm =( . u m (A7) 

V Vim J \ IVVlm J 

It is easy to check that B tm and Ci m satisfy the following orthogonality relations 

(B tj ,C lm )=0, (A8a) 

(Bij, Bi m ) = (Cij, Clm) = aim • $il • Sjm, (A8b) 

in which 5ij is Kronecker-5. The velocity components can now be expressed in terms of 
vector spherical harmonics as 
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N I ' 

a) = i Vim W m cos mX + ^ m sin mX ) 

1=0 m=0 

+ Wi m {-c\ m cos mX + 4 m sin m\) } . (A10) 



The coefficients, 6 r m , fef 71 , c," 1 , c- m , can be evaluated using the following inverse trans- 
forms: 



b\ m = — <| / dcj) cos cj)Vi m I dXvcosmX — I d(f> cos 4>Wi m I dXusmmX 

y/Otlm I J-TL JO J-TL Jo 



Jm 



(All) 



1 J rf y^r yf /-27T 

6 ?; m = ——= <( / c/0 cos <p Wi m I dXu cos mX + / cos VJ m / <iA v sin mX 



o ./o 



(A12) 



c( is <-/> V nm J dXu cos mA + J d(f> cos Wz m J dX v sin mA 



0413) 



Jm 



COS </> Wi 



I in 



277 



(iA v cos mX + 



cos V h 



2- 



dX u sin mX 



(AU) 



Substituting (A9-A14) in e quations (7a) and (7b ) and performing a few steps of algebra 
(see equations A. 14 and A. 15 of ISwarztrauber I (119961 ) for detailed derivation), we can obtain 
the following expressions for the variables ( and 5, which no longer contain the unbounded 
terms that cause the pole problem in spherical polar coordinate system: 
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Similarly following ISwarztrauber I ( 119961 ) (see his equations A. 16 and A. 17) the expression 
for the gradient of the scalar variable, h, can be given as 
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By grouping the derivative and algebraic terms and using the compact notations, defined 



as 



E u = (C + 2w c sin0)f, 
E v = — (£ + 2u; c sin0)M, 
the time-evolution equations for u and t> can be written as 
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Note that E u and E v can be expressed in terms of V\ m and W\ m in an analogous way as 
u and v have been expressed in the Equations A9 and A10. 

In order to derive the time-evolution equations for u and v in spectral space we differ- 
entiate equations (A11-A14) with respect to time and obtain the following: 
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in which, u and v respectively denote the time derivatives of u and v. Substituting u and v 
from equations (A21) and (A22) in the above equations (A23-A26), we obtain 
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cos0(9A \ 2 



+ / d0 cos V/ m / rfA sin mA • — ( h G/i 

7-5 Jo d( t> 
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b\ m (s) 



AA, AW t AX X 1 9 (U 2 + V\ 

d<p cos0 Wi m \ dX cosmA • 7 -^- ( — h Gn 



cos 4> dX 



+ / a0 cos0 Vz m / dX smmX ■ — | h Gn 

7-; Jo 



(A28) 



dc', 



lm 



dt i/ofr. 



cos 6 V n 



r-2-K 

Jo 



dX E u cosmX + 



2tt 



cos Wi m I dX E v sin mX 



f 2lT 1 d (u 2 + v 2 

dip cos0 Vi m I dX cosmA ■ -— ( h Gh 

J cos (p oX 



f 2lT d (u 2 + v 2 
+ / d(p cos W^ m / dA sin mX • 77- ( — h G/i 



<90 



= 4 m (s) 



1 



i-2-k 

/ ' 

JO 



dcj) cos V/ m / dX cos mA 



1 d fu 2 + v 



+ dS cos W^ m / dX sin mA • — 
lo d<t> 



2n 



cos <9A 
<9 f u 2 + v 2 



+ Gh 



+ Gh 



(A29) 



g c lm 



p2tt 



p2w 

Jo ' 



d(j) cos Wj m / dX E v cos mA + / d(f> cos VJ m / dX E u sin mA 

7T 

9 f U 2 + V 2 



y/<Xlr, 



/ ' 

Jo 



o?0 cos W^ m / o?A cos mA • — 

ocp 



+ Gh 



— J d(f) cos Vi m / <iA sin mA 
J -S Jo 



1 d f u 2 + v 



c (s) - 



v 7 ^ 



cos <9A 



2tt 



cos Wj m / dA cos mA • — — 

O0 



+ G7i 



9 [ u 2 + v 2 



+ Gh 



p2n 



d<{) cos VJ m / <iA sin mA 



1 d fu 2 + v 



cos 9 A 



+ Gh 



(A30) 



As in the case of the meteorological primitive equations, the shallow water equations in 
the solar case admit high frequency gravity waves as well as low frequency shearing flows. 
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In order to increase the efficiency of the computation, we employ an implicit scheme for the 
linear terms and an explicit scheme for the nonlinear terms of the Equations (A27-A30). 
Therefore we sort out those terms and represent them in their corresponding spectral forms. 



Starting with the Equation (A27), we have the following linear term: 



/ dd cos (f)Vi m / dXcosmX G— — / deb cos dWi m I dXsinmX G 
J-Z Jo V 9(f) J J_* J V 



1 dh 

cos 6 dX 



Substituting the expressions for ^^f| and || from Equations (A17) and (A18) into the 
above expression, and applying the orthogonality relation from Equation (A7), we obtain 



1 



/2tt ____ 
dXcosmX ^2 (h + ^)GV hmi (a l r imi cos mi A + a- 1 " 11 sin mi A) 
li=0mi=0 



N h 



- I 2 d<j> cos <f>W lm ( rfAsinmA^ J2 V(h + l )GW hmi ■ (a[ lfni cos mi A - a l r imi sin mi A) 

"'"f ^° (i=0mi=0 

iV !i * .1 

= V V 4 im V^ (ii + i) / d<fi cos • c (v^ imi + W, ro Wi imi ) 



-27T 



!i=0mi=0 



-a l ™Gy/i (/ + !). 



(A31) 



A part of the nonlinear terms of equation (A27) contains the derivative of (u 2 + v 2 )/2. 
Noting that (u 2 + v 2 )/2 is a scalar term which is represented in spectral space as 



u 2 + 



TV I 



= ^2 ^2 Pim ( dl ™ cos m ^ + ^ m sin ' 



(A32) 



Z=0 m=0 



we can use similar formulae as given in expressions (A17) and (A18) for the derivative of a 
scalar variable to evaluate the derivatives of (u 2 + v 2 )/2. After a few steps of algebra, those 
terms can be simplified to 

/-27T Q / 2 _|_ y 2 

d(f) cos 4> Vim I dX cos mA • — 



/f 
rf0 COS W im 



2 " M ,19 /« 2 + t; 2 

ctA sin mA • 



o 



cos 4> dX 
= -(G-a l r m + d l r m ) (/ + !). 
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Performing similar analysis for the nonlinear terms involving derivatives of {u 2 + v 2 )/2 
respectively in equations (A28-A30), we get 



f 2n 1 d fu 2 + v 2 

dip cos (f)Wi m / dXcosmX- - — 

J Q COS CM 



f 2lT d (u 2 + v 2 

I d(f> cos Vi m J ^ sinm ^'7^( — 2 — 

= -{g- a[ m + d[ m ) y/Hl + Tj, 



cos <p Vim I dX cos mX 







cos <9A 



r2ir d ( u 2 + 

+ / dtp cos Wj m / dX sin mA • — 

" 2 



0, 



1-2-K 

/ 1 



cl0 cos Wj m / c/A cos mX ■ — 



d ( u 2 + 



rf r 2 - 
— / dcj) cos VJ m / dX sin mA 
./-i Jo 



90 V 2 
1 d fu 2 + v 2 



cos <9A 



0. 



Recalling Equation (3) from §2.1, we now write it as 

dh 



dt 



in which 



u dh dh 
$ = _ 



(A33), 



(A34) 



COS CM 90 

To obtain the time-evolution equations for the spectral coefficients involved in h, we differ- 
entiate Equations (A3) and (A4) with respect to time and using equation (A33) we obtain 



dai m 



dt ^/af r . 



d<t>Pi 



tm 



i>2tt 
(0)COS0 / 

Jo 



dX cosmA • \P 



v 7 ^ 



1 r n 
/ d(j)Pi m (0) cos / dX cos mX ■ 5, 

0Li m J_* J 



(A35) 
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da[ m 1 



d<pPi m ((f)) cos 4> d\ sin mX ■ \1> (0, A) 
Jo 



dt y/afr, 

1 f 2lT 
/ d<f>Pi m ((f)) cos / <iAsinmA-<5. (A36) 

\JOLlm J-TL Jo 

Substituting in equations (A35) and (A36), the expression for 8 from equation (A16), 
and using orthogonality of Legendre polynomials we get 



dal m 1 



dt y/afr 



[ 2 d<f>Pi m ((f)) cos (j> f dX cosmA • \P + y/l (I + l)b l ™, 
J-l Jo 

= a l r m (V) + (I + l)b l r m , (A37) 
d(f)P lm ((f)) cos / dX sin mX ■ # (0, A) + ^/Z (Z + l)^ m - 

a( m (^) + v 7 ^ (I + l)6| m , (^38) 



da im in 



The Equations (A27-A30) and (A37-A38) represent the full set of hydrodynamic shallow- 
water Equations (1-3) of §2.1 in the spectral form. Like the Equations (A37) and (A38), the 
Equations (A27-A30) can also be expressed in compact form as 



= b l r m (S) - ^/l(l + l)d l ; n - Gy/l (l + l)a l r m , (A39) 

db lm 

= b\ m (S) - ^HT+Y)d l :r - Gy/T(lTl)a<r, (A40) 

c l r m (S) . (AA1) 



dt 



Note that due to certain symmetry properties (antisymmetry here) of Vi m and Wi m , 2nd 
and 3rd terms in the right hand side of equation (Equation (29)) get cancelled. Similar 
cancellation happens also in Equation (30) for which is reproduced in compact form 

below, namely 

^ = c\ m (S) . (AA2) 
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In order to evaluate the integrals efficiently and accurately, we use Gaussian grids in both 
latitude and longitude directions. Denoting the non- uniformly spaced Gaussian grid points 
in (j) by 4>j an d the corresponding Gaussian weights by Wj, a ^-integral can be expressed as 

7T 

/ (0) = f 2 d<Pf (fa « Yl w if M ■ ( A43 ) 
Similarly a A-integral respectively involving a cosine and sine functions can be expressed as: 

A (A) = / dX cos mXf (X) &'S^ w k cos mX k f (X k ) , (A44) 
Jo 



C 

h(X)= / d\ sinmXf (A) « w fc sinmA fe / (A fc ) , (A45) 

in which A& is the Gaussian grid in longitude and w k is the corresponding weighting function. 
Using the expressions (A43-A45), the expansion coefficients for the height and the flow 
variables can be represented in Gaussian grids as: 



a 1 ™ = i WjPi m ((f)j) cos fa w k cos mX k h A fc ) , (AA6) 

v j k=l 

l N k 

a[ m = ——= V] WjPim cos fa y~]w k sinmX k h (fa, X k ) , (AA7) 

b l r m = Wj COS fa { Vi m ((f) j ) VCm (<t>j) - W lm ((f) j ) US m ((f) j) } , (A48) 

j 

b[ m = W 3 COS <i>j i W lm (fa) UC m ((f) j) + Vi m ((f) j) VS m (fa)}i, (A49) 

j 

c l r m = Y ™j cos fa {V lm ((f) j) uc m ((f) j) + W lm ((f) j) vs m ((f) j)} , (A50) 
j 

^ = Yl W i C ° S ^3 {- W lm (fa) VC m (fa) + V lm ((f) j) us m (fa)} , (A51) 



in which 



uc m 

(<A?) = ^2 cos (^fc m ) M (fa, A fc ) , 
k=i 

us m (fa) = Y sin(X k m)u (fa, X k ) , 



k=i 
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vc 



vs 



k=l 



k=l 



Swarztrauber I (119961 ) solved the shallow-water equations in physical space evaluating 



nonlinear terms in spe ctral space. We will be instead evolve the spectral coefficients following 
Hack fc Jakob I (119921 ). so we construct a semi- implicit time evolution scheme, because the 
time step utilized by explicit time differencing methods is limited by the speed of the fastest 
gravity wave. The semi- implicit time integration removes this constraint by treating terms 
that give rise to gravity waves implicitly, and the remaining terms explicitly. 

Each of the Equations (A37-A42) has a linear term. We write all of these equations in 
centrally time differenced form and the linear term averaged over time r — 1 and r + 1. We 
omit the superscript of lm from all the terms for focussing only on time evolution scheme. 



2A T 



[K +1 



br 1 ] 



b T r (E) - ^n{n + l)d T r - G^n{n + 1) 



< +1 + a T r 



T-l 



(A52) 



2A T 



[bl +1 - br 1 ] = bl(E) - y/n(n + 1R - Gy/n{n + l) 



< +1 + aJ~ 1 



(A53) 



2A, 



[a T r +1 - a;- 1 } = + y/n(n + l) 



Vr+l 



+ 61 



T-l 



2A 



(A54) 



(A55) 



Equations (A52) and (A54) are coupled, similarly Equations (A53) and (A55) are also 
coupled. We write Equations (A52) and (A54) in matrix form as 



gA T ^n (n + 1) 1 

1 -A rV /n (n+T) 



2A T 61 (E) 



n 




(n + l)d T r - gA T ^n (n + l)^ -1 + 6 



r 



2A T a T r (*) + A T y/n (n + l)^" 1 + ap 1 



(A56) 
(A57) 
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which can be written as 

1 



a r - = T - XJ - 7 — TTV — n [a; 1 {-gA 2 T n{n+l) + l} 



[gA 2 T n (n + l) + 1] 

+2A r - A r n(n + 1)< + (fop 1 + A T 6^(S)) Vn(n + 1)}] (A58) 

-2A r {-K(E) + « + 9 <~' + gAXW) x/n(n+l)}] (A59) 
Similarly, Equations (A53) and (A55) can be written as 

< +1 = b^^+D + n [<-'{- g AM» + D + i} 

+2A r {<(*) - A T n(n + + (ft,""- 1 + A r 6[(3)) V"(n + 1)}] (A60) 

-2A r {-6J(H) + (d* +jAr<(*)) + (A61) 

We discretize the Equations (A41) and (A42) in time for explicit evolution: 

c T r +1 = c T r + A^p 1 (S) (A62) 
= cI + A^rMS) (A63) 

We have presented above all the equations we need to do time integration. At each 
time step we integrate equations (A58-63). The details of the time evolution in each time- 
advancement involve sequentially the following steps: 

1. From 6j. m , 6- m , cj, m and c' m compute u and v in (6>, 0) Gaussian grid using Equations 
(A9) and (A10). 

2. From a 1 ™ and a- m compute ft, in (9,<j)) Gaussian grid using Equation (Al). 

3. From b l ™, b\ m , c 1 ™ and cf 11 compute ( and 5 in (0, 0) Gaussian grid using Equations 
(A15) and (A16). 

4. From u, v, ( and 5 computed in steps 1 and 2, and u c given, compute E u and E v using 
Equations (A19) and (A20). 

5. Substitute E u and S„, obtained from step 4, in Equations (A11-A14) with -u and t> 
respectively to evaluate b l ™(E), &- m (H), 4 m (S) and c- m (S). 
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6. Compute |(w 2 + v 2 ) from the already computed values of u and v from step 1, and 
then using sequentially the Equations (A32), (A17), (A18), (A3) and (A4), compute 
d l r m and df l . 

7. Compute the gradients of \(u 2 + v 2 ) in (#,</>) Gaussian grids using d l r m and <i- m from 
step 6, using Equations (A17) and (A18). 

8. Compute the gradients of h in (9, 0) Gaussian grids using aj™ and a 1 ™ in Equations 
(A17) and (A18). 

The above steps constitute the complete computations of right hand side of Equations 
(A37-A42). Using the right hand side we compute b 1 ™, a 1 ™, 6- m , a- m for the next time 
advancement from Equations (A58-A61), and cj. m and c- m from Equations (A62) and (A63). 
We repeat the above eight steps forward as long as we want to continue the time marching. 
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